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Abstract 

Random dense packing of poly-disperse discs and spheres generated under the 
infinitesimal gravity protocol is investigated by a Monte Carlo simulation. For 
the binary discs packed in two dimensions, the packing fraction of disc assembly 
becomes lower than that of the monodisperse system when the size ratio is close 
to unity. We show that the suppressed packing fraction is caused by an increase 
of the adjacent neighbours with long bonds where the adjacent neighbours is 
defined on the basis of the Laguerre (radical) tessellation. 

For the poly-disperse systems in two and three dimensions, the packing frac- 
tion is shown to have a minimuma as a function of the poly-dispersity. Percola- 
tion process in the densely packed discs and spheres is also studied. The critical 
area (volume) fraction in two (three) dimensions is shown to be a monotonically 
increasing (decreasing) function of the poly-dispersity. 

1 Introduction 

Random close packing of spheres is an important concept in understanding the struc- 
ture and properties of liquids and glasses. The first significant contribution in this 
direction was reported by BernalpQ who analyzed the structure of glasses on the basis 
of the sphere packing. It is also well known that solidification of liquids [2 j and the flow 
of granular materials 3j are closely related to packing properties of hard spheres. 

In practical applications, it is common to encounter systems in which the size of 
dispersed particles is not unique but has a certain distribution. For example, a func- 
tional composite system can be fabricated by depositing nano-clusters which are syn- 
thesized by the plasma-gas-condensation technique. In this case, the size distribution 
of the nano-clusters plays important role in determining the properties of the composite 
system[4j. Beside the packing fraction, the connectivity of particles is also important 
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characteristics of the system composed with dispersed spheres. For example, a print- 
ing ink consists of carbon blacks dispersed in a varnish, and the network formation of 
carbon blacks is one of the key parameters controlling the properties of the ink. It is, 
therefore, important in practical applications to find out the dependence of packing 
fraction and connectivity on the width of the distribution or the poly-dispersity. 

In this paper, we investigate the packing fraction for a random dense packing of bi- 
nary hard discs and of poly-disperse hard discs and spheres. In order to avoid subtlety 
in the definition of random close packing[5], we study the dense packing produced by 
the drop-and-roll or the infinitesimal gravity protocol using a Monte Carlo simulation. 
Although the structure produced by this protocol is not the random close packing nor 
the maximally random jammed stated EI] , it is worth studying the structure produced 
by this sequential protocol since in practical applications the packing of particles is 
produced by a fixed procedure such as molecular vapour deposition or dispersing parti- 
cles in a solvent. In §2 we explain briefly the protocol employed. The packing of binary 
discs in two dimensions is studied in §3, where results of computer simulation and an 
analysis based on the Laguerre tessellation are given. In §4, we present results for the 
packing of poly-disperse discs in two dimensions and spheres in three dimensions. For 
the poly-disperse systems we also study the percolation process of randomly selected 
discs and spheres in the packed structure. Conclusion is given in §5. 

2 Infinitesimal gravity protocol 

In the infinitesimal gravity protocoljT, a container of L x L (L x L x L) is prepared 
in two (three) dimensional space. A particle (a disc in two dimensions and a sphere in 
three dimensions) is introduced far above the container at a horizontal position selected 
randomly. The radius of the particle is chosen so that the distribution of the radius 
R obeys a given function (j)(R). Then, the particle is dropped vertically towards the 
bottom of the container. The particle drops freely or rolls down around a particle that 
it touches until the particle settles into a stable position or makes contact with the 
base line. (See Fig.l.) We repeat the process until the container is filled with particles. 
It is known that this protocol produces for a mono-disperse system a structure very 
close to the random close packing, though the structure is not the maximally jammed 
oneO Ej. 

3 Binary discs 

3.1 Monte Carlo simulation 

Using the infinitesimal gravity protocol, we produced packed structures of two kinds 
of discs, one with radius Rs and the other with radius Rl- The system size is set 
to L/2Rl = 100 and discs at given fraction x of the smaller disc are packed up to a 
height of 140 x 2Rl. To reduce the effects of the boundaries, upper and lower 20 x 2Rl 
layers are discarded (so that the remaining area is 100 x 100 x 4i?|), and we imposed 
a periodic boundary condition in the horizontal direction. Then, the packing fraction 
/ is obtained as a function of the disc size ratio r = Rs / Rl and the area fraction of 
smaller discs n = xR 2 s /[{l — x)R\ + xR\\. 
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Figure 1: The infinitesimal gravity protocol is illustrated schematically. After contact- 
ing with another disc, each disc, A and B, rolls down around the target disc maintaining 
contact until it makes contact with another disc. It rolls down around the new disc if 
it can (A), and if it cannot roll down further, then it settles there (B). 

The packing fractions f(n,r) obtained from the Monte Carlo simulation are illus- 
trated in Fig. 2. Note that /(0,r) = /(l,r) = f(n, 1) = /(n,0) = f 2 ^ 0.82. Here, 
f% is the packing fraction of the random close packing of monodisperse discs [8 a . For 
r ~ 0, as n is increased from n = 0, f(n,r) increases dramatically, because of the va- 
cant spaces between larger discs being filled by smaller discs. The increase of packing 
fraction reaches a maximum at a certain n where the vacant spaces are almost filled up, 
and decreases again to the value for the monodisperse assembly. On the other hand, 
f(n, r) for r ~ 1 decreases from / 2 as n is increased from n = 0, and is almost constant 
over the wide range of n until it increases again to f 2 near n = 1. The decrease of 
f(n,r) is due to the fact that the smaller disc, which is too large to fill the vacant 
spaces, replaces the larger discs and destroys the structure formed by larger discs [7j. 

It is interesting to compare the contour for the packing fraction with the contour 
of the constant mean squared deviation of radii scaled by the mean diameter, ((-R 2 ) — 
(i?) 2 )/4(i?) 2 , in the (n, r) plane, which is shown in Fig. 3. Here, (R) and (R 2 ) are the 
mean and the mean squared radius of discs in the system, respectively. The similarity 
of these contours suggests that the packing fraction is determined mostly by the relative 
width of the distribution function. This observation is strongly supported by the results 
for poly-disperse systems discussed in §4. 

3.2 Analysis by the Laguerre tessellation 

Using the Laguerre tessellation 0, we can define adjacent neighbours of a disc and 
connect centres of all adjacent neighbour pairs to form a triangular network which fills 
the space without gaps. 

For an infinite system, the mean number of adjacent neighbours for each disc is 
equal to 6, which is a consequence of Euler's condition, and therefore the number of 
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Figure 2: A 3D plot of the packing fraction / of a binary hard disc assembly as 
a function of size ratio r and area fraction n. Each value was determined from an 
average of 20 samples. The contours are drawn from 0.81 to 0.85 every 0.05. 



triangles included in the system is twice as many as the number of discs. As a result, 
the packing fraction of the system / is expressed in term of the mean area of triangles 

A m as 

f = a) 

We analysed the bond lengths in detail. The bond length distribution in our simu- 
lation is shown in Fig. 4 for the bond between two adjacentlarge discs. It is apparent 
from Fig. 4 that the major difference of bond distributions is in the behaviour of the 
long bond region. We can conclude that, although the tail is small compared with the 
whole distribution, the behaviour of the tail greatly affects the packing fraction[7j. 



4 Poly-disperse discs and spheres 
4.1 Packing fraction 

We examined the Gaussian and the uniform distributions for two and three dimensions. 
For two dimensions, we also examined the log-normal distribution. All distribution 
functions are normalized so that the average radius J R<f)(R)dR is equal to m and 
the mean squared deviation f(R — m) 2 (p( y R)dR is equal to a 2 . We set the average 
diameter 2m as the unit of length and use the width a as a parameter representing the 
poly-dispersity. 

We prepared a 200 x 200 (60 x 60 x 60) square (cubic) box and examined 1000 
samples in two dimensions and 3000 samples in three dimensions to determined the 
packing fraction. In this simulation, fixed boundary conditions are imposed in hori- 
zontal directions. 

Figures 5 and 6 show the dependence of / on the poly-dispersity a. The packing 
fraction / = 0.8175 ± 0.0002 at a = in two dimensions is consistent with the value 
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Figure 3: Contour plot of the constant mean squared deviation of radii of discs. Con- 
tours are shown for logo" = —2 ~ for every 0.1, where a = J (R 2 ) — (R) 2 /2(R) . 



f 2 = 0.82 reported in literature jHj- Similarly, / = 0.5575 ± 0.0003 at a = for three 
dimensions agrees well with the value fs = 0.555 ± 0.005 reported by Onoda et a/.[TU]. 
When a is increased, the packing fraction / decreases and the dependence is roughly 
independent of the distribution function. When a is increased further, the packing 
fraction / exhibits a minimum and increases again. This behaviour is identical to the 
behaviour of the packing fraction of binary discs discussed in §3. 

4.2 Continuum percolation 

In order to construct percolation process in the random assembly of poly-disperse 
discs and spheres, we randomly select a given fraction of the total number of par- 
ticles (coloured white at the beginning), irrespective of their size, and change their 
colour to black. Two black particles are regarded as connected when they touch each 
other. Then, we investigate the percolation of the black particles when their fraction 
is increased. 

We determined the critical percolation area (volume) fraction A c (V c ) of black par- 
ticles as a function of a, where percolation was judged by the appearence of a cluster 
spanning from the bottom to the top surface. In Figs. 7 and 8, A c and V c are plotted as 
a function of a. Our value A c = 0.4777 ± 0.0005 and V c = 0.1938 ± 0.0001 at a = are 
close to the critical area fraction determined from lattice models A c = 0.45 ± 0.02 
and the value of V c = 0.185 ± 0.005 in literature [HJ. 

The critical area fraction A c is an increasing function of a and is roughly inde- 
pendent of the distribution function when a is small. In fact, if we use a power-law 
function 

A c {a) = a 2 + b 2 a c * (2) 
to fit the data for 0^<t^0.15, parameters are give by 

a 2 = 0.4777 ± 0.0005, b 2 = 0.028 ± 0.002, c 2 = 0.54 ± 0.05. (3) 

The solid line in Fig. 7 represents Eq. ((2)). 
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Figure 4: The distribution function of bond length between two large discs at r = 0.9. 
The solid, dash-dot and dot curves correspond ton = 0.0 (the mono-diserse assembly), 
n — 0.1 and n = 0.5. Note that each distribution function has a delta function-like 
peak at the length which equals 2Rl. 

The critical volume fraction V c is a decreasing function of a. Although the de- 
pendence is weak, this is an amazing contrast to the critical area fraction A c in two 
dimensions. For 0<<r<0.15, V c can well be fitted by the following line 

V c (a) =a 3 + b 3 a (4) 

with 

a 3 = 0.1938 ±0.0001, b 3 = -0.005 ± 0.001, (5) 

which is shown by the solid line in Fig. 8. 

The significant difference in the poly-dispersity dependence of the critical area and 
volume fractions can be explained by the surface area of discs and spheres which 
determines the effectiveness of connection P.3J. In fact, we can easily show that while 
the average circumference of discs is independent of a, the average surface area of 
spheres increases as a is increased. 

5 Concluding remarks 

We have examined by Monte Carlo simulation the packing fraction of poly-disperse 
discs and spheres which are packed in two and three dimensions by the infinitesimal 
gravity protocol. It is a striking result that the packing fraction decreases for the weak 
poly-dispersity. This is due to the fact that a little increase of poly-dispersity makes 
particles hard to adjust vacant spaces and further increase of a introduces smaller 
particles filling up vacant spaces. It is an open and interesting mathematical question 
if one can rigorously prove this effect in general. 

We have aslo studied the conitinuum percolation of poly-disperse systems in two and 
three dimensions and have shown that the critical area fraction in two dimensions and 




Figure 5: The relation between / and a for various size distribution functions in two 
dimensions: The results are shown for Gaussian (circles), uniform (squares), and log- 
normal (triangles) distributions. 

the critical volume fraction in three dimensions have opposit dependence on the poly- 
dispersity. Furthermore, the critical area (volume) fraction is shown not to depend 
much on the distribution function. These results indicate that the concept of the 
critical area (volume) fraction as a dimensional invariant has to be applied with 
some caution. 

This work was supported in part by the grant from the Ministry of Education, 
Science, Sports and Culture. 
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Figure 6: The relation between / and a for two size distribution functions in three 
dimensions: The results are shown for Gaussian (circles) and uniform (squares) distri- 
butions. 
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Figure 7: The critical area fraction A c is plotted against a for three size distribution 
functions: Gaussian (circles), uniform (squares) and log-normal (triangles). The solid 
line represents the fitting function (J2J). Error bars are much smaller than the size of 
the symbols. 
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Figure 8: The critical volume fraction V c is plotted against a for two size distribution 
functions: Gaussian (circles) and uniform (squares). The solid line represents function 
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